Ideas discussed at last meeting:
First look at fitting models for different numbers of cells. With only oocyte and nurse cell 2 we get:
source('run_model_comparison_stst.R')
out = run_model_comparison_stst(identifier='Ststv203',use_real_data=TRUE,run_mcmc=FALSE,nSamples=6,nTest=14,parametersToPlot = c('nu','xi','phi'),verbose=FALSE,show_diagnostic_plots=FALSE,test_on_mutant_data=FALSE, cell_indices = c(1,2,rep(0,14)))
## Loading required package: ggplot2
## Loading required package: StanHeaders
## rstan (Version 2.17.3, GitRev: 2e1f913d3ca3)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
## Loading required package: tidyr
##
## Attaching package: 'tidyr'
## The following object is masked from 'package:rstan':
##
## extract
## Saving 7 x 5 in image
## Warning: Removed 154 rows containing missing values (geom_point).
With only oocyte and neighbouring nurse cells (1,2,3,5,9):
out = run_model_comparison_stst(identifier='Ststv204',use_real_data=TRUE,run_mcmc=FALSE,nSamples=6,nTest=14,parametersToPlot = c('nu','xi','phi'),verbose=FALSE,show_diagnostic_plots=FALSE,test_on_mutant_data=FALSE, cell_indices = c(1,2,3,5,9,rep(0,11)))
## Saving 7 x 5 in image
## Warning: Removed 154 rows containing missing values (geom_point).
With only oocyte and nurse cells within two teps from the oocyte (1,2,3,4,5,7,8,9,10,11,13):
out = run_model_comparison_stst(identifier='Ststv205',use_real_data=TRUE,run_mcmc=FALSE,nSamples=6,nTest=14,parametersToPlot = c('nu','xi','phi'),verbose=FALSE,show_diagnostic_plots=FALSE,test_on_mutant_data=FALSE, cell_indices = c(1,2,3,4,5,7,8,9,10,11,13,rep(0,5)))
## Saving 7 x 5 in image
## Warning: Removed 154 rows containing missing values (geom_point).
And with all the nurse cells in the tissue:
out = run_model_comparison_stst(identifier='Ststv206',use_real_data=TRUE,run_mcmc=FALSE,nSamples=6,nTest=14,parametersToPlot = c('nu','xi','phi'),verbose=FALSE,show_diagnostic_plots=FALSE,test_on_mutant_data=FALSE, cell_indices = seq_len(16))
## Saving 7 x 5 in image
## Warning: Removed 154 rows containing missing values (geom_point).
Just with the cells further away from the oocyte:
out = run_model_comparison_stst(identifier='Ststv207',use_real_data=TRUE,run_mcmc=FALSE,nSamples=6,nTest=14,parametersToPlot = c('nu','xi','phi'),verbose=FALSE,show_diagnostic_plots=FALSE,test_on_mutant_data=FALSE, cell_indices = c(rep(0,11),6,12,14,15,16))
## Saving 7 x 5 in image
## Warning: Removed 154 rows containing missing values (geom_point).
Everything except the neighbouring nurse cells:
out = run_model_comparison_stst(identifier='Ststv208',use_real_data=TRUE,run_mcmc=FALSE,nSamples=6,nTest=14,parametersToPlot = c('nu','xi','phi'),verbose=FALSE,show_diagnostic_plots=FALSE,test_on_mutant_data=FALSE, cell_indices = c(rep(0,5),4,6,7,8,10,11,12,13,14,15,16))
## Saving 7 x 5 in image
## Warning: Removed 154 rows containing missing values (geom_point).
Now see if we can summarise the distribution of values for \(\nu\) that we get for fitting to each different choice of cells
library(dplyr)
library(ggplot2)
stst_df = data.frame(nu=numeric(),xi=numeric(),phi=numeric(),which_cells=integer())
for (id in seq(from=3,to=8)){
print(paste('fits/alt_save_MCStstv20',id,'.Rsave',sep=''))
load(paste('fits/alt_save_MCStstv20',id,'.Rsave',sep='')) #load in the results of fitting at stst for each choice of cells
e$which_cells = rep(id,length(e$nu))
stst_df = full_join(stst_df,as.data.frame(e))
}
## [1] "fits/alt_save_MCStstv203.Rsave"
## Joining, by = c("nu", "xi", "phi", "which_cells")
## Warning: Column `nu` has different attributes on LHS and RHS of join
## Warning: Column `xi` has different attributes on LHS and RHS of join
## Warning: Column `phi` has different attributes on LHS and RHS of join
## [1] "fits/alt_save_MCStstv204.Rsave"
## Joining, by = c("nu", "xi", "phi", "which_cells")
## Warning: Column `nu` has different attributes on LHS and RHS of join
## Warning: Column `xi` has different attributes on LHS and RHS of join
## Warning: Column `phi` has different attributes on LHS and RHS of join
## [1] "fits/alt_save_MCStstv205.Rsave"
## Joining, by = c("nu", "xi", "phi", "which_cells")
## Warning: Column `nu` has different attributes on LHS and RHS of join
## Warning: Column `xi` has different attributes on LHS and RHS of join
## Warning: Column `phi` has different attributes on LHS and RHS of join
## [1] "fits/alt_save_MCStstv206.Rsave"
## Joining, by = c("nu", "xi", "phi", "which_cells")
## Warning: Column `nu` has different attributes on LHS and RHS of join
## Warning: Column `xi` has different attributes on LHS and RHS of join
## Warning: Column `phi` has different attributes on LHS and RHS of join
## [1] "fits/alt_save_MCStstv207.Rsave"
## Joining, by = c("nu", "xi", "phi", "which_cells")
## Warning: Column `nu` has different attributes on LHS and RHS of join
## Warning: Column `xi` has different attributes on LHS and RHS of join
## Warning: Column `phi` has different attributes on LHS and RHS of join
## [1] "fits/alt_save_MCStstv208.Rsave"
## Joining, by = c("nu", "xi", "phi", "which_cells")
## Warning: Column `nu` has different attributes on LHS and RHS of join
## Warning: Column `xi` has different attributes on LHS and RHS of join
## Warning: Column `phi` has different attributes on LHS and RHS of join
stst_df = stst_df %>% mutate(which_cells=factor(which_cells))
h <- ggplot(data = stst_df,aes(x=nu,color=which_cells)) + geom_density() + theme_bw()
h
This could certainly be interpreted as evidence of non-uniform bias in transport through ring canals spatially across the tissue.
** Two super-compartments
To reduce stochastic variation due to small sample sizes, we further consider agreggating into two super-compartments. I assume well mixed concentrations of mRNA in each super-compartment:
\[\begin{array} \frac{\text{d}x_1} = 7a - b*(1-\nu)x_1 + b\nu x_2, \\ \frac{\text{d}x_1} = 8a - b\nu x_2 + b(1-\nu) x_1. \end{array}\]From the long time limit, we get \[ x_2 = \frac{8a + b(1-\nu) x_1}{b\nu}. \]
library(rstan)
library(mvtnorm)
library(dplyr)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
#read in data on WT
N=6 # sample size
data = matrix(as.numeric(read.csv('data/exp_data.csv',sep=',',header=FALSE,stringsAsFactors = FALSE)),ncol=16,byrow=TRUE)
data = data[seq_len(N),]
#aggregate into two compartments
x1_obs = data[,seq(from=1,to=16,by=2)] %>% apply(.,1,sum)
x2_obs = data[,seq(from=2,to=16,by=2)] %>% apply(.,1,sum)
estimates <- stan(file = 'super_compartment_model.stan',
data = list (
N = N,
x1_obs = x1_obs,
x2_obs = x2_obs
),
seed = 42,
chains = 4,
warmup = 1000,
iter = 2000,
control = list(adapt_delta = 0.99)
)
print(estimates, pars = c('a','b','theta','sigma'))
## Inference for Stan model: super_compartment_model.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## a 34.62 0.74 26.65 1.50 13.67 28.63 49.36 101.01 1288 1
## b 79.57 1.32 61.30 3.79 31.63 64.86 114.72 225.94 2145 1
## theta 0.44 0.00 0.12 0.27 0.35 0.41 0.49 0.74 1241 1
## sigma 247.92 1.13 44.43 173.60 215.69 244.27 274.91 348.34 1547 1
##
## Samples were drawn using NUTS(diag_e) at Thu Aug 16 13:30:06 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).
pairs(estimates, pars = c('a','b','theta','sigma'))
The results from the two super-compartment model contradict what we have found using the full data suggesting that transport through ring canals should be less biased than previously thought. However, these result rely upon an assumption of well mixed concentrations of mRNA within each super-compartment. This assumption does not hold. I would suggest that the results here may be misleading due to the inaccuracy of this assumption.
TODO: try forward simulating with different models, but particularly with different \(\nu\) for single connection to oocyte. Plus make predictions for WT and OE from full fits.
Another possibility: assume the oocyte is an absorbing state. Currently we assume bidirectional transport through each ring canal. We could hypothesize that the connections directly to the oocyte are `special’ in some way such that transport only occurs in one direction here. This should decrease the concentration of mRNA in the neighbouring nurse cells and increase in the oocyte in the model.
Might be plausible for the full dynamic model.
But at steady state, the absorbing oocyte case results in qualitatively different behaviour in the sense that for this case the eigenvector corresponding to the largest eigenvalue is [1,0 … 0] irrespective of the transport bias \(\nu\). mRNA accumulates in the oocyte but not in the nurse cells, whereas assuming all the ring canals are bidirectional results in a distribution of mRNA across cells dependent on \(\nu\) with nonzero values in the nurse cells.
library(rstan)
mc <- stan_model('model_comparison5.stan')
expose_stan_functions(mc) #get hold of functions defined in the stan code
nu = 0.72 #pick a vlaue for the bias from ring canals
B = construct_matrix(0.72)
set_oocyte_absorbing = function(B){
B[,1] = 0
return(B)
}
B = construct_matrix(nu)
B_absorb = set_oocyte_absorbing(B)
print(eigen(B))
## eigen() decomposition
## $values
## [1] -1.917756e+00 -1.663930e+00 -1.493222e+00 -1.374166e+00 -1.202139e+00
## [6] -1.163084e+00 -1.101435e+00 -1.000000e+00 -7.461741e-01 -7.102357e-01
## [11] -6.698615e-01 -6.258343e-01 -5.080965e-01 -4.627503e-01 -3.613152e-01
## [16] 1.976128e-16
##
## $vectors
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.666340488 -0.663418548 -0.53797776 -0.664653592 0.46128567
## [2,] -0.666340488 -0.107998368 0.53797776 -0.108199422 -0.46128567
## [3,] -0.181435015 0.663418548 -0.32163923 -0.108199422 -0.12560151
## [4,] 0.181435015 0.107998368 0.32163923 -0.017613859 0.12560151
## [5,] -0.116039538 0.180639412 0.24850174 0.664653592 0.43121093
## [6,] 0.116039538 0.029406416 -0.24850174 0.108199422 -0.43121093
## [7,] 0.031595912 -0.180639412 0.14857102 0.108199422 -0.11741259
## [8,] -0.031595912 -0.029406416 -0.14857102 0.017613859 0.11741259
## [9,] -0.088193947 0.115530699 0.11894222 0.180975697 -0.18875355
## [10,] 0.088193947 0.018807323 -0.11894222 0.029461160 0.18875355
## [11,] 0.024013954 -0.115530699 0.07111164 0.029461160 0.05139490
## [12,] -0.024013954 -0.018807323 -0.07111164 0.004796003 -0.05139490
## [13,] 0.015358492 -0.031457362 -0.05494158 -0.180975697 -0.17644726
## [14,] -0.015358492 -0.005120966 0.05494158 -0.029461160 0.17644726
## [15,] -0.004181898 0.031457362 -0.03284776 -0.029461160 0.04804407
## [16,] 0.004181898 0.005120966 0.03284776 -0.004796003 -0.04804407
## [,6] [,7] [,8] [,9] [,10]
## [1,] 0.54381084 0.49517497 0.679901590 -0.42402411 0.49193414
## [2,] 0.08852735 -0.49517497 0.110681654 0.42402411 0.08008230
## [3,] -0.54381084 0.29604884 0.110681654 0.11545572 -0.49193414
## [4,] -0.08852735 -0.29604884 0.018017944 -0.11545572 -0.08008230
## [5,] 0.32512663 0.17451114 0.110681654 0.07384147 -0.13394665
## [6,] 0.05292759 -0.17451114 0.018017944 -0.07384147 -0.02180527
## [7,] -0.32512663 0.10433448 0.018017944 -0.02010598 0.13394665
## [8,] -0.05292759 -0.10433448 0.002933154 0.02010598 0.02180527
## [9,] -0.25119615 -0.28713514 -0.679901590 -0.52152792 0.45986119
## [10,] -0.04089240 0.28713514 -0.110681654 0.52152792 0.07486112
## [11,] 0.25119615 -0.17166866 -0.110681654 0.14200462 -0.45986119
## [12,] 0.04089240 0.17166866 -0.018017944 -0.14200462 -0.07486112
## [13,] -0.15018192 -0.10119308 -0.110681654 0.09082122 -0.12521365
## [14,] -0.02444822 0.10119308 -0.018017944 -0.09082122 -0.02038362
## [15,] 0.15018192 -0.06050002 -0.018017944 -0.02472932 0.12521365
## [16,] 0.02444822 0.06050002 -0.002933154 0.02472932 0.02038362
## [,11] [,12] [,13] [,14] [,15]
## [1,] 0.44367490 0.591241031 -0.46310664 0.56496562 0.55109740
## [2,] -0.44367490 0.096248540 0.46310664 0.09197115 -0.55109740
## [3,] 0.26525864 0.096248540 0.12609734 -0.56496562 0.32948302
## [4,] -0.26525864 0.015668367 -0.12609734 -0.09197115 -0.32948302
## [5,] -0.20494152 -0.591241031 -0.43291317 0.33777438 0.19421950
## [6,] 0.20494152 -0.096248540 0.43291317 0.05498653 -0.19421950
## [7,] -0.12252780 -0.096248540 0.11787609 -0.33777438 0.11611746
## [8,] 0.12252780 -0.015668367 -0.11787609 -0.05498653 -0.11611746
## [9,] 0.32668014 0.353483583 -0.18424064 0.19910698 0.15471424
## [10,] -0.32668014 0.057543839 0.18424064 0.03241276 -0.15471424
## [11,] 0.19531132 0.057543839 0.05016610 -0.19910698 0.09249856
## [12,] -0.19531132 0.009367602 -0.05016610 -0.03241276 -0.09249856
## [13,] -0.15089951 -0.353483583 -0.17222858 0.11903952 0.05452489
## [14,] 0.15089951 -0.057543839 0.17222858 0.01937853 -0.05452489
## [15,] -0.09021786 -0.057543839 0.04689539 -0.11903952 0.03259864
## [16,] 0.09021786 -0.009367602 -0.04689539 -0.01937853 -0.03259864
## [,16]
## [1,] -0.9490332012
## [2,] -0.1544937769
## [3,] -0.1544937769
## [4,] -0.0251501497
## [5,] -0.1544937769
## [6,] -0.0251501497
## [7,] -0.0251501497
## [8,] -0.0040942104
## [9,] -0.1544937769
## [10,] -0.0251501497
## [11,] -0.0251501497
## [12,] -0.0040942104
## [13,] -0.0251501497
## [14,] -0.0040942104
## [15,] -0.0040942104
## [16,] -0.0006664994
print(eigen(B_absorb))
## eigen() decomposition
## $values
## [1] -1.8506469 -1.5891102 -1.4773990 -1.2839774 -1.1946481 -1.1423873
## [7] -1.0921857 -0.8600000 -0.7384633 -0.6948672 -0.6649883 -0.5760226
## [13] -0.5002191 -0.4336353 -0.3414496 0.0000000
##
## $vectors
## [,1] [,2] [,3] [,4]
## [1,] 3.943199e-01 -4.505475e-01 4.150432e-01 5.366668e-01
## [2,] -8.485429e-01 -4.571711e-15 -7.130051e-01 1.305241e-15
## [3,] 3.227778e-16 8.325228e-01 -3.233952e-16 -1.398704e-16
## [4,] 2.801942e-01 -1.001964e-15 -3.515094e-01 -3.824201e-16
## [5,] -6.917153e-18 6.434661e-17 1.423595e-16 -8.012420e-01
## [6,] 1.629329e-01 1.387343e-15 3.534887e-01 7.874473e-16
## [7,] -5.719377e-17 -2.749042e-01 -4.177991e-17 -7.839515e-16
## [8,] -5.380145e-02 1.134009e-15 1.742688e-01 -5.957913e-16
## [9,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [10,] 1.199176e-01 1.340689e-15 1.616794e-01 -4.123511e-16
## [11,] -4.485057e-17 -1.598567e-01 6.855239e-17 9.528992e-16
## [12,] -3.959755e-02 3.382269e-17 7.970747e-02 3.808785e-16
## [13,] 8.247258e-19 -1.243206e-17 -2.817746e-17 2.645751e-01
## [14,] -2.302596e-02 2.603537e-16 -8.015629e-02 -4.154894e-16
## [15,] 7.879598e-18 5.278570e-02 1.595270e-17 4.849430e-16
## [16,] 7.603318e-03 -9.846383e-16 -3.951681e-02 2.591564e-16
## [,5] [,6] [,7] [,8] [,9]
## [1,] 4.334712e-01 5.176061e-01 -4.982278e-01 -0.7071068 -5.799736e-01
## [2,] -6.021460e-01 -1.377579e-14 6.327410e-01 0.0000000 4.980107e-01
## [3,] 2.960621e-15 -6.875658e-01 -3.279721e-16 0.0000000 2.987860e-16
## [4,] 1.988324e-01 -7.902666e-16 3.119394e-01 0.0000000 -1.644462e-01
## [5,] 1.782552e-16 1.590424e-16 3.467737e-18 0.0000000 -9.722521e-18
## [6,] -5.105009e-01 -7.471665e-15 2.077652e-01 0.0000000 -9.562545e-02
## [7,] 1.187878e-15 -3.389679e-01 -2.080447e-16 0.0000000 -3.342942e-17
## [8,] 1.685706e-01 1.748434e-15 1.024276e-01 0.0000000 3.157613e-02
## [9,] 0.000000e+00 0.000000e+00 0.000000e+00 0.7071068 0.000000e+00
## [10,] 2.519077e-01 5.936144e-15 -3.815211e-01 0.0000000 5.736660e-01
## [11,] -1.240031e-15 3.408765e-01 2.133050e-16 0.0000000 2.289866e-16
## [12,] -8.318151e-02 4.790898e-16 -1.880887e-01 0.0000000 -1.894281e-01
## [13,] -9.758457e-17 -9.307494e-18 -4.009275e-18 0.0000000 1.332493e-17
## [14,] 2.135680e-01 3.309307e-15 -1.252753e-01 0.0000000 -1.101524e-01
## [15,] -5.446099e-16 1.680511e-01 1.113667e-16 0.0000000 8.644795e-17
## [16,] -7.052150e-02 -1.089499e-16 -6.176032e-02 0.0000000 3.637301e-02
## [,10] [,11] [,12] [,13]
## [1,] 6.674936e-01 6.450812e-01 -8.012420e-01 7.575452e-01
## [2,] -8.803343e-15 -4.988041e-01 1.398727e-16 -4.406262e-01
## [3,] -5.393249e-01 -3.306966e-17 6.086083e-16 1.434219e-15
## [4,] 1.320414e-15 -2.459089e-01 6.343025e-16 1.454976e-01
## [5,] -6.896561e-17 -1.493518e-16 5.366668e-01 -3.212961e-16
## [6,] -7.532866e-16 2.472936e-01 -5.930919e-16 -3.735641e-01
## [7,] 1.780885e-01 -2.218286e-16 -5.989054e-16 4.644482e-16
## [8,] 2.087871e-15 1.219150e-01 -8.141743e-17 1.233532e-01
## [9,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [10,] -6.653132e-15 -3.580943e-01 7.763897e-17 -1.714590e-01
## [11,] -4.572410e-01 7.301365e-17 9.532322e-16 5.286705e-16
## [12,] 5.623556e-17 -1.765394e-01 4.809293e-16 5.661684e-02
## [13,] 4.604350e-17 -1.450831e-16 2.645751e-01 -1.067148e-16
## [14,] 1.040675e-15 1.775335e-01 -3.130082e-16 -1.453634e-01
## [15,] 1.509839e-01 -2.285894e-16 -3.722299e-16 2.178551e-16
## [16,] 7.177216e-16 8.752347e-02 -1.399036e-16 4.799990e-02
## [,14] [,15] [,16]
## [1,] 8.606270e-01 -9.006089e-01 1
## [2,] 6.511794e-15 3.575727e-01 0
## [3,] -4.339514e-01 -3.215760e-16 0
## [4,] -3.316949e-16 1.762823e-01 0
## [5,] -2.695983e-16 -1.068387e-17 0
## [6,] 4.292391e-15 1.174116e-01 0
## [7,] -2.139367e-01 -1.348413e-16 0
## [8,] -9.261474e-16 5.788358e-02 0
## [9,] 0.000000e+00 0.000000e+00 0
## [10,] 2.344028e-15 9.653869e-02 0
## [11,] -1.424911e-01 -8.822899e-17 0
## [12,] -4.093421e-16 4.759328e-02 0
## [13,] -6.069709e-17 -1.802243e-18 0
## [14,] 1.762836e-15 3.169919e-02 0
## [15,] -7.024770e-02 -3.379160e-17 0
## [16,] -5.409596e-16 1.562760e-02 0
Now try to visualise the effect of changing \(\nu\) by forward simulating from the model. (For this we will assume bidirectional transport through all ring canals again).
library(tidyr)
#keep all other parameters constant during this experiment
m0 = c(0, rep(1,15)) #initial condition
th = c(4,300) #c(96.5,249)
sig = 50
phi = 0.289
nSamples = 15
nTest = 5
nTotal = nSamples + nTest
source('extract_times_and_scaling.R')
times = extract_times_and_scaling(nSamples,nTest,test_on_mutant_data=FALSE)
data = matrix(as.numeric(read.csv('data/exp_data.csv',sep=',',header=FALSE,stringsAsFactors = FALSE)),ncol=16,byrow=TRUE)
raw_data = data[times$sort_indices2,]
#take a vector of transport bias values for nu to loop through
nu_vec = seq(from=0.7, to=1, by=0.05)
all_extracted_samples = data.frame(rna=numeric(),time=numeric())
for (j in seq_along(nu_vec)){
B = construct_matrix(nu_vec[j])
samples <- stan(file = 'mrna_transport6.stan',
data = list (
T = nTotal,
y0 = m0,
t0 = times$t0,
ts = times$ts2,
theta = array(th, dim = 2),
sigma = sig,
phi = phi,
B = B %>% as.vector()
),
algorithm="Fixed_param",
seed = 42,
chains = 1,
iter =100,
refresh = -1
)
estimates = rstan::extract(samples, pars='y_hat')
#make the samples go in a 'nice' format
pred <- as.data.frame(estimates, pars = 'y_hat') %>%
gather(factor_key = TRUE) %>%
group_by(key) %>%
summarize(lb = quantile(value, probs = 0.05),
median = quantile(value, probs = 0.5),
ub = quantile(value, probs = 0.95))
xdata <- data.frame(rna = as.vector(raw_data),cellID = as.vector(matrix(rep(1:16,nTotal),nrow=nTotal,byrow=TRUE)),time = rep(times$ts2,16))
pred <- pred %>% bind_cols(xdata) %>%
mutate(split = if_else(time %in% times$ts3,'train','test'),
nu = nu_vec[j])
all_extracted_samples = full_join(all_extracted_samples, pred)
}
##
## Elapsed Time: 0 seconds (Warm-up)
## 0.04661 seconds (Sampling)
## 0.04661 seconds (Total)
## Joining, by = c("rna", "time")
##
## Elapsed Time: 0 seconds (Warm-up)
## 0.04929 seconds (Sampling)
## 0.04929 seconds (Total)
## Joining, by = c("rna", "time", "key", "lb", "median", "ub", "cellID", "split", "nu")
##
## Elapsed Time: 0 seconds (Warm-up)
## 0.052461 seconds (Sampling)
## 0.052461 seconds (Total)
## Joining, by = c("rna", "time", "key", "lb", "median", "ub", "cellID", "split", "nu")
##
## Elapsed Time: 0 seconds (Warm-up)
## 0.056709 seconds (Sampling)
## 0.056709 seconds (Total)
## Joining, by = c("rna", "time", "key", "lb", "median", "ub", "cellID", "split", "nu")
##
## Elapsed Time: 0 seconds (Warm-up)
## 0.054649 seconds (Sampling)
## 0.054649 seconds (Total)
## Joining, by = c("rna", "time", "key", "lb", "median", "ub", "cellID", "split", "nu")
##
## Elapsed Time: 0 seconds (Warm-up)
## 0.059533 seconds (Sampling)
## 0.059533 seconds (Total)
## Joining, by = c("rna", "time", "key", "lb", "median", "ub", "cellID", "split", "nu")
##
## Elapsed Time: 0 seconds (Warm-up)
## 0.055752 seconds (Sampling)
## 0.055752 seconds (Total)
## Joining, by = c("rna", "time", "key", "lb", "median", "ub", "cellID", "split", "nu")
p1 <- ggplot(all_extracted_samples, aes(x = time, y = median, group = factor(nu), color=factor(nu)))
p1 <- p1 + geom_line() +
facet_wrap(~cellID,scales='free') + #needed to remove factor(cellID)
labs(x = "time", y = "rna") +
theme(text = element_text(size = 12), axis.text = element_text(size = 12),
strip.text = element_text(size = 8)) +
scale_color_discrete(name = "nu") +
scale_colour_brewer(palette=7) +
geom_point(aes(x = time, y = rna))
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
p1
## Warning: Removed 1078 rows containing missing values (geom_point).
And finally compare this to what happens when we change one the transport bias for ring canals to the oocyte, with bias \(\nu\) for the other ring canals fixed at \(0.9\).
B1 = construct_matrix(0.72) #we fix ring canals for all other nurse cells
#take a vector of transport bias values for nu to loop through
nu_vec = seq(from=0.7, to=1, by=0.05)
all_extracted_samples = data.frame(rna=numeric(),time=numeric())
for (j in seq_along(nu_vec)){
#mess around with ring canals into oocyte
B2 = construct_matrix(nu_vec[j])
B = B1
B[,1] = B2[,1]
B[1,] = B2[1,]
diag(B) = diag(B) - apply(B,2,sum)
samples <- stan(file = 'mrna_transport6.stan',
data = list (
T = nTotal,
y0 = m0,
t0 = times$t0,
ts = times$ts2,
theta = array(th, dim = 2),
sigma = sig,
phi = phi,
B = B %>% as.vector()
),
algorithm="Fixed_param",
seed = 42,
chains = 1,
iter =100,
refresh = -1
)
estimates = rstan::extract(samples, pars='y_hat')
#make the samples go in a 'nice' format
pred <- as.data.frame(estimates, pars = 'y_hat') %>%
gather(factor_key = TRUE) %>%
group_by(key) %>%
summarize(lb = quantile(value, probs = 0.05),
median = quantile(value, probs = 0.5),
ub = quantile(value, probs = 0.95))
xdata <- data.frame(rna = as.vector(raw_data),cellID = as.vector(matrix(rep(1:16,nTotal),nrow=nTotal,byrow=TRUE)),time = rep(times$ts2,16))
pred <- pred %>% bind_cols(xdata) %>%
mutate(split = if_else(time %in% times$ts3,'train','test'),
nu = nu_vec[j])
all_extracted_samples = full_join(all_extracted_samples, pred)
}
##
## Elapsed Time: 0 seconds (Warm-up)
## 0.057958 seconds (Sampling)
## 0.057958 seconds (Total)
## Joining, by = c("rna", "time")
##
## Elapsed Time: 0 seconds (Warm-up)
## 0.052128 seconds (Sampling)
## 0.052128 seconds (Total)
## Joining, by = c("rna", "time", "key", "lb", "median", "ub", "cellID", "split", "nu")
##
## Elapsed Time: 0 seconds (Warm-up)
## 0.0512 seconds (Sampling)
## 0.0512 seconds (Total)
## Joining, by = c("rna", "time", "key", "lb", "median", "ub", "cellID", "split", "nu")
##
## Elapsed Time: 0 seconds (Warm-up)
## 0.051745 seconds (Sampling)
## 0.051745 seconds (Total)
## Joining, by = c("rna", "time", "key", "lb", "median", "ub", "cellID", "split", "nu")
##
## Elapsed Time: 0 seconds (Warm-up)
## 0.052531 seconds (Sampling)
## 0.052531 seconds (Total)
## Joining, by = c("rna", "time", "key", "lb", "median", "ub", "cellID", "split", "nu")
##
## Elapsed Time: 0 seconds (Warm-up)
## 0.049678 seconds (Sampling)
## 0.049678 seconds (Total)
## Joining, by = c("rna", "time", "key", "lb", "median", "ub", "cellID", "split", "nu")
##
## Elapsed Time: 0 seconds (Warm-up)
## 0.049038 seconds (Sampling)
## 0.049038 seconds (Total)
## Joining, by = c("rna", "time", "key", "lb", "median", "ub", "cellID", "split", "nu")
p2 <- ggplot(all_extracted_samples, aes(x = time, y = median, group = factor(nu), color=factor(nu)))
p2 <- p2 + geom_line() +
facet_wrap(~cellID,scales='free') + #needed to remove factor(cellID)
labs(x = "time", y = "rna") +
theme(text = element_text(size = 12), axis.text = element_text(size = 12),
strip.text = element_text(size = 8)) +
scale_color_discrete(name = "nu") +
scale_colour_brewer(palette=7) +
geom_point(aes(x = time, y = rna))
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.
p2
## Warning: Removed 1078 rows containing missing values (geom_point).
Changing the transport bias, \(\nu\), for selected ring canals neighbouring the oocyte has limited effect relative to predictions from the simpler model where \(\nu\) is assumed to be spatially uniform.
Seems possible though to look at fitting this model with spatially varying \(\nu\) values based on whether cells are neighbouring the oocyte. Can then compare between the models to assess whether this is plausible.
source('mrna_transport_full.R')
mrna_transport_inference_full('full_spatial_nu_v201', use_real_data = TRUE, run_mcmc = FALSE, nSamples=6, nTest=14, show_diagnostic_plots = TRUE, test_on_mutant_data = FALSE)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## [1,] 566 13 18 NA 16 NA NA NA 10 NA NA NA NA
## [2,] 708 47 117 50 91 34 86 49 76 94 100 63 110
## [3,] 800 44 119 NA 29 NA NA NA 5 NA NA NA NA
## [4,] 1106 162 36 118 92 131 61 76 82 129 39 135 122
## [5,] 446 83 22 NA 2 NA NA NA 44 NA NA NA NA
## [6,] 982 77 69 NA 15 NA NA NA 0 NA NA NA NA
## [7,] 1185 75 122 NA 51 NA NA NA 23 NA NA NA NA
## [8,] 1076 22 19 NA 67 NA NA NA 70 NA NA NA NA
## [9,] 1007 51 5 NA 32 NA NA NA 11 NA NA NA NA
## [10,] 1032 31 144 NA 43 NA NA NA 6 NA NA NA NA
## [11,] 1078 47 159 NA 18 NA NA NA 61 NA NA NA NA
## [12,] 1454 261 144 180 108 142 10 96 18 127 20 53 19
## [13,] 1012 23 141 NA 28 NA NA NA 88 NA NA NA NA
## [14,] 1862 48 32 NA 73 NA NA NA 21 NA NA NA NA
## [15,] 1045 78 130 NA 101 NA NA NA 11 NA NA NA NA
## [16,] 1840 283 194 168 190 179 67 22 85 70 56 47 174
## [17,] 1365 102 318 NA 155 NA NA NA 31 NA NA NA NA
## [18,] 1615 98 80 NA 284 NA NA NA 164 NA NA NA NA
## [19,] 1485 407 53 384 399 174 47 152 21 315 37 176 266
## [20,] 1652 68 376 249 43 166 28 330 263 293 208 73 29
## [,14] [,15] [,16]
## [1,] NA NA NA
## [2,] 73 81 71
## [3,] NA NA NA
## [4,] 84 31 35
## [5,] NA NA NA
## [6,] NA NA NA
## [7,] NA NA NA
## [8,] NA NA NA
## [9,] NA NA NA
## [10,] NA NA NA
## [11,] NA NA NA
## [12,] 109 20 11
## [13,] NA NA NA
## [14,] NA NA NA
## [15,] NA NA NA
## [16,] 136 20 17
## [17,] NA NA NA
## [18,] NA NA NA
## [19,] 40 99 73
## [20,] 272 12 98
## Warning: Removed 154 rows containing missing values (geom_point).
## Saving 7 x 5 in image
## Warning: Removed 154 rows containing missing values (geom_point).
## Joining, by = "parameter"
## Saving 7 x 5 in image
## This is bayesplot version 1.6.0
## - Online documentation and vignettes at mc-stan.org/bayesplot
## - bayesplot theme set to bayesplot::theme_default()
## * Does _not_ affect other ggplot2 plots
## * See ?bayesplot_theme_set for details on theme setting
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## the following parameters were dropped because they are duplicative
## a b
## Warning: Removed 154 rows containing missing values (geom_point).
Finally, make predictions for full model for WT and OE, with spatially uniform \(\nu\) and spatially varying \(\nu\).